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Abstract. Recently Baumann et al. [arXiv:0709.3228vl] studied the phase-ordering 
kinetics of the two-dimensional Ising model for T < T c with uniform spatially quenched 
disorder by Monte-Carlo simulations. They found that the two-time response and 
correlation functions are in agreement with the predictions of local scale invariance 
C3 ■ generalised to Z 7^ 2. The present paper shows why this is actually not true and 

suggests an alternative approach which leads to a much better agreement with the 
numerical results. 

O 

£\\ , Submitted to: Journal of Statistical Mechanics: Theory and Experiment 

> 

PACS numbers: 05.50.+q, 05.70.Ln, 64.60.Ht 

(N 

<N 

1. Introduction 

<^) A ferromagnetic systems starting with a disordered initial state exhibits phase ordering 

kinetics, forming ordered domains whose linear size grows with time as £ ~ t 1 '*, where 
z is the dynamical exponent. For homogeneous systems with a non-conserved order 
parameter one has z = 2 [1], while in systems with spatially quenched disorder the 
dynamical exponent is found to vary continuously with the intensity of the disorder. 
The probably simplest example of such a system is the two-dimensional ferromagnetic 
Ising model with quenched disorder [2]. It is defined by the Hamiltonian 

H = — 2_^ JijO~iO-j , ctj = ±1 , (1) 

where the sum runs over all nearest neighbours of a two-dimensional square lattice with 
periodic boundary conditions. The coupling constants J^ = J^ are uncorrelated random 
variables distributed uniformly over the interval 

J„£[l-e/2, l+e/2], (2) 



£ 



Response function of the disordered Ising model 2 

where e G [0, 2] is a parameter controlling the strength of the disorder. The model evolves 
by conventional heat bath dynamics, i.e., for each local update a site i is randomly 
selected and the spin o~i is oriented according to the local field hi = Y^uj) Jij a j °f the 
neighboring spins by setting 

{h ■ IhT 
+1 with probability eh . /fc r^ e _ h . /fcT ^ 

—1 otherwise. 

The model exhibits a second-order phase transition between a paramagnetic and a 
ferromagnetic phase at some critical temperature T c (e) > which depends on the 
intensity of the disorder. In the ferromagnetic phase T < T c (e), starting with a random 
initial configuration, one observed phase ordering with a dynamical exponent z > 2. 
According to phenomeno logical scaling arguments [2] and field-theoretic studies [3,4] its 
value is given bjtjj 

^ = 2 + ^. (4) 

For example, for T = 1, - a temperature well below T c that will be considered throughout 
this paper - the dynamical exponent varies between 2 and 4 depending on the strength 
of the disorder e. 

Phase ordering in disordered ferromagnets gives rise to physical ageing, charac- 
terised by a non-exponential relaxation, dynamical scaling, and breaking of time- 
translational invariance [6]. An interesting quantity, which allows one to study ageing 
phenomena, is the response function 



R(t,S^: 'W^) 



5h(s,0) 



(5) 



which describes the averaged spatio-temporal response of the magnetisation to a local 
perturbation by a weak magnetic field h. The response function can be measured by 
applying a weak field at a particular site at time s and monitoring subsequent spin 
flips caused by the field at distance r at some later time t > s. Since in disordered 
ferromagnetic Ising model with heat bath dynamics a positive external field h(s, 0) > 
can only turn spins up but never down, the response function is strictly positive in the 
present case. 

Ageing phenomena are characterised by dynamical scaling. This postulate implies 
that the response function obeys the scaling form 

HC^O-.-'-wg,^;), (6) 

| Note that this relation is still debated in the literature. In particular, Huse and Henley [5] conjectured 
that domains grow as a power of logt. However, the simulations performed in the present work seem 
to be consistent with Eq. 0. 
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where 1Z is a scaling function and 

a = l/z. (7) 

During the past decade there have been various attempts to generalise the concept of 
dynamical and anisotropic scaling to some kind of local scale invariance [7], in a similar 
sense as conformal invariance generalises global scale invariance of two-dimensional 
equilibrium phase transitions to local scale transformations. Such theories of local scale 
invariance (LSI) together with explicit representations of its generators make non-trivial 
predictions about the specific form of scaling functions. Recently, LSI was developed 
further as to include systems with z ^ 2 [8] . As a first application to a non-linear model 
with z t^ 2, Baumann et al. [9] applied this generalised version of LSI to the disordered 
ferromagnetic Ising model, finding the predictions in good agreement with the numerical 
results. In contrast to these results, the aim of the present work is 

(a) to show that the predictions of generalised LSI, as quoted in [9], are in principle 
incompatible with the actual response function of the disordered Ising model, 

(b) to demonstrate this incompatibility numerically by a direct measurement of the 
response function, 

(c) to suggest an alternative approach with algebraically distributed waiting times 
which yields consistent results in much better agreement with the numerical data. 

The present results for T < T c complement previous studies Pleimling, Gambassi, 
Corberi, Lippiello and Zannetti [10, 11] who have questioned the applicability of LSI 
to the critical Ising model without quenched disorder at criticality. A similar discussion 
concerning the application of LSI to critical directed percolation without quenched 
disorder can be found in Refs. [12, 13]. 

2. Test of the predictions of LSI 

According to Henkel et al. [8,9,14] the generalised version of LSI makes the following 
predictions: 

(i) The scaling function TZ factorises into two parts 

<? (i^vH* GM(^) • (8) 

where we may choose .F(O) = 1. 

(ii) The first factor in this expression, which determines the form of the autoresponse 
R(t, s) = s~ 1 ~ a fii(t/s), again factorises into two power laws of the form 
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where r is a proportionality constant, a' is another exponent which is generally 
different from a, and Xr is the so-called autoresponse exponent. 

(iii) The second factor, which determines the spatial profile of the response, can be 
expressed as an integral 

^ P \u) = J-0^\kfeM^-k-a\k\ z ) (10) 

parametrised by a dimensioned non-universal parameter a and another exponent j3. 

Baumann et al. tested these predictions indirectly by measuring the integrated 
spatio-temporal response and autocorrelation function, finding a convincing agreement. 
However, this result is surprising because the integral reminds of a Levy-stable 
distribution generated by anomalous diffusion. Such distributions are expected to occur 
in systems with algebraically distributed long-range interactions in space, typically in 
the range 1 < z < 2. The present model, however, does not fall into this category 
because all interactions are local. It is this disturbing discrepancy which motivated the 
present work. 

In the following we demonstrate that integral expression ( fTOl) is in principle incom- 
patible with the response function of the disordered Ising model in 2D evolving by heat 
bath dynamics. To see this, one has to evaluate the integral ffTUl) in two dimensions 
using polar coordinates 

J* a >®(&) = -L / dkkk?e- akZ i d^e iuhcos4. 



o Jo 

oo 



' ' dkk 1+f} e- akz J (uk), (11) 



JT(a,(3) { 



a * 



256a 



U 



4vr j 

where Jo denotes the Bessel function of the first kind. The remaining integral can be 
evaluated exactly for certain rational values of z. As an example let us consider the 
special case z — 4. Here we obtain the result 

-^(^^3 (1 + 1; 1,1,1; 
16vr 
„»a-|- 1 r(| + l) 1 F»(S + l;l,f,f; 2r)I 

647 " • ^ 

where p F q ({ai, . . . , a p } ; {b±, . . . , b q } ; z) are generalised hypergeometric functions. Inser- 
ting the numerically determined values a = 0.15(2) and (3 = —0.04(10) reported in [9], 
this function displays an oscillatory behaviour, becoming negative for certain values of it. 
However, as mentioned before, the response function of the disordered ferromagnetic 
Ising model evolving by heat bath dynamics is strictly positive. Hence we can conclude 
that the predictions of LSI, as expressed by Eq. ffTUl) . cannot be applied to the disordered 
Ising model studied here. 
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Figure 1. Prediction of LSI for the scaling function ^^'^'(u) for a = 0.15 
and (3 = —0.04 (blue line) compared with numerical data of simulations at 
temperature T = 1 below T c (red dots) for z = 4. The numerical profile has 
been obtained by using the algorithm described in Appendix A, with s = 512, 
t = 32768 on a system with N = 750 2 sites. 



This incompatibility can be supported by a direct numerical measurement of the 
response profile. As already observed in the homogeneous kinetic Ising model [10, 11, 
15,16], the full response function resolves the deviations from LSI much better than the 
integrated response (see Appendices for technical details). Fig. [1] shows the numerically 
measured profile and the analytical profile predicted by LSI. As can be seen, the profiles 
differ significantly. 

3. An alternative approach 



In the following I suggest an alternative approach which is guided by phenomenological 
arguments in the spirit of Paul, Puri and Rieger [2]. Starting point is the observation 
that it makes a big difference whether a local perturbation by an external field flips 
a spin in the interior of an ordered domain or a spin adjacent to a domain wall. In 
the interior of a domain the generated minority island shrinks extremely fast driven 
by surface tension so that its contribution to the response function is extremely short- 
ranged in time and can be neglected. Spin flips at the boundaries, however, may displace 
a domain wall and thus have a much longer life time. In other words, the response to 
a local perturbation is always confined to the domain walls and thus propagates in the 
same way as the domain walls move in space. Without quenched disorder, the domain 
walls propagate and merge diffusively, hence the local perturbation spreads in space like 
a random walk, leading to a Gaussian profile. 

Quenched disorder changes the dynamics of the domain walls drastically. In 
particular, domain walls in a disordered environment may be pinned along lines of 
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strongly coupled sites [2]. Such a pinned domain wall remains immobile for a certain 
waiting time. In a scale-free system it is reasonable to assume that these waiting times 
At are algebraically distributed as -P(At) ~ (At) _1_K with k < 1, leading effectively 
to a subdiffusive dynamics. With this physical picture in mind it is near at hand to 
conjecture that the response to a local perturbation in a disordered ferromagnet spreads 
essentially in the same way as a short-ranged random walk with algebraically distributed 
waiting times between subsequent random moves. 

To solve this problem, let us first introduce a parameter r which is proportional 
to the number of diffusive moves performed so far. Initially the process starts with 
t = t = 0. As time proceeds one may ask the question how the elapsed time t is 
distributed for given r. This distribution P{t\r) evolves by according to the equation 

(d« + d T )p{t\T) = (13) 

with the initial condition P{t\0) = o~(t), where d£ denotes a temporal fractional 
derivative which is known to generate algebraically distributed waiting times (see 
e.g. [17,18] and references therein). To solve this equation consider the Fourier transform 



P(t\r) = -;= / du e iuJt P(lu\t) (14) 

by which Eq. ([TBI is turned into 

(^) k + 9 t )p(cj|t) = (15) 



with the solution P{oj\t) = e T ( tuJ ) K . By means of the inverse Fourier transform one 
arrives at 

1 f°° 
P(t\r) = —= / du e iwt ~ r{iu)) \ (16) 

where we identify k = 2/z by dimensional analysis. This integral can be solved only in 
certain special cases, e.g. 



P(t\r) 



' S(t-r) ifz = 2 

■<«.MM^ + ! MtHa: i[z = 3 (17) 



2v / 37ri 7 / 3 2 v / 37rf 5 / 

e ~^t 
2y^i 3 / 2 



if 2 = 4 



Note that this solution is already properly normalised by J °° dt P(t\r) = 1 and obeys 
the required initial condition P(t\0) = S(t). 

The solution gives us the distribution of the physical time t after a given number 
of diffusive moves r. However, to compute the response profile one needs to know 
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how many diffusive moves have been performed at a given time t. This conditional 
distribution P(r\t) is related to P(t\r) by 

t 



P(r\t) = ±^dt>P(t>\T 



giving 



P{r\t) 



' S(t-r) if = 2 

if z = 3 (19) 



4V37Tt 4 /3 

e ^T 



if 2 = 4 



We now assume that deep in the ageing regime <> s the normalised response is just 
a superposition of Gaussians with different widths according to the varying number of 
diffusive moves performed so far, i.e. 



R(t,s,r) 
R(t,s,0) 



POO 

/ dre- r2/DT P(r\t). (20) 

Jo 



Here D is the diffusion constant of the random walker which here plays the role of a 
non-universal metric factor. Unlike the LSI prediction discussed in the previous section, 
this function is strictly positive. For the special case z = 4 this integral can be evaluated, 
giving 

2^-y3,0 / r 4 I Inn 

R(t,s,r) r G 0,3 ^iStl -2' U ' U , 

R{t,s,0) K A^DVt ' 

where G'T'f 1 denotes the Meijer G-function. This function can be evaluated numerically 
and may be compared with the normalised numerical data, where D plays the role of a 
fit parameter. As shown in Fig. [2], both curves are in convincing agreement for D = 2.9. 

In order to demonstrate that the obtained profile differs from a Gaussian, the 
integral (1201) may be approximated by the saddle-point method. This yields the 
expression 

- 4/3 ' ' v^/3 



R(t,s,0) V3 { ' 

which decays far away from the origin as 

, R(t,s,r) x 4 / 3 

- ln W^oj ~ W (23) 

confirming non-Gaussian behaviour. As shown in the inset of Fig. [21 the exponent 4/3 
is in fair agreement with the numerical results. 
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Figure 2. Profile according to Eq. (|2ip (blue) compared to the same numerical 
data as in Fig. [I] (red) . The inset illustrates the asymptotic decay, confirming 
the power 4/3 predicted in Eq. (|23|) . 



4. Summary 



The work by Baumann et al. was guided by the expectation that LSI generalised to z ^ 2 
should be a generic theory for ageing phenomena. The disordered Ising model seemed 
to be a natural candidate to test this expectation in the range 2 < z < 4. Studying 
the integrated response, the authors found an almost perfect coincidence. However, as 
demonstrated in the present work, a direct measurement of the response function reveals 
significant deviations. Moreover, the predicted response function becomes negative for 
certain distances while the actual response of the Ising model has to be strictly positive. 

In the second part of the paper I demonstrated that another fractional differential 
equation, which is local in space but non-local in time, gives results in fair agreement 
with the numerical data. I do not claim that this differential equation has any 
fundamental meaning, the only purpose is to demonstrate that a simple approach guided 
by physical intuition in the spirit of Ref. [2] leads to a much better agreement with the 
numerical data. Another possible explanation in terms of super-ageing was proposed 
by Paul et al. in Ref. [19] and it would be interesting to investigate how closely these 
approaches are related. The question whether it is possible to construct a different 
representation of LSI describing the disordered Ising model remains open. 

Finally it should be mentioned the present findings do not support a recently 
suggested hypothesis of super- universality [20], stating that domain growth in 
homogeneous and weakly disordered Ising systems is described by identical scaling forms. 

Acknowledgements: 
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Appendix A. Detailed description of the simulation algorithm 

Previous studies on the homogeneous Ising model have shown that the integrated 
response is not a suitable quantity to detect deviations from LSI, instead one has to 
measure the full response function. A direct measurement of the response function is 
notoriously difficult to perform because of strong fluctuations. As a way out, Chatelain 
proposed to relate the response function to a certain correlation function measured in a 
system without external field [21]. Based on this idea several authors devised powerful 
algorithms to measure the response function in the Ising model [11,16,22]. 

Here I present an alternative algorithm for the Ising model with standard heat 
bath dynamics, where perturbations at all sites are processed simultaneously - a trick 
which allows one to obtain reasonable statistical averages within short time. A similar 
algorithm was already used for the contact process (see Ref. [12]). In the present case the 
algorithm is somewhat more complicated, as will be explained in the following. However, 
the algorithm is conceptually simpler as previous approaches in so far as it measures 
the response function directly without relying on the non-trivial proof by Chatelain. 

The standard kinetic Ising model with heat bath dynamics is defined as follows. 
Each site % of a d- dimensional lattice is associated with classical variables o~i — ±1 which 
represents the orientation of a local spin. Initially the system is prepared in a disordered 
state, where all spins independently take the values ±1 with equal probability. Then, 
starting at time t — 0, the system evolves at finite temperature (3 = 1/ksT by random- 
sequential updates according to the following dynamical rules: 

• Select a random site i. 

• Calculate the local field hi = ^ • o~j, where j runs over the 2d nearest neighbours 
of site i, respecting the chosen boundary conditions. 

• Generate a random number z G (0, 1) and set 

(new) ._ f +1 if z < |(l + tanh(/3/ij)) 
I — 1 otherwise 

As usual in models with random sequential updates, the time variable t is incremented 
by 1/N, where N = L d denotes the total number of sites. 

Obviously, the third step can be reformulated equivalently as follows: 

• Define a threshold field 

h\ c) ■= /5~ 1 atanh(22 - 1) , (A.2) 

where z G (0, 1) is a random number, and set 

a (new) ;= j +1 \ih i >hf , A3 , 

1—1 otherwise 
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As we will see below, this formulation of the update rule is more convenient for a 
measurement of the autoresponse function. 

In order to measure the autoresponse function R(t, s) directly we would have to 
apply a weak external field < h\ f x <^C 1 at time s locally at a randomly chosen site i 
and to study how much the ensemble average of the local magnetisation (o~j(t)) at a 
different site j changes at some later time t > s. However, such a measurement is hard 
to perform since a localised weak field leads only to occasional spin flips in the future so 
that the average response to such a single spin flip is hard to quantify. This is probably 
the reason why most researches avoid a direct measurement of the autoresponse function 
and prefer to study the integrated response. In the kinetic Ising model, however, there 
are two ways to optimise a direct measurement, namely by (i) reweighting the probability 
of spin flips and (ii) simultaneous tagging of all spins in the system. 

The algorithm suggested here seems to have a comparable performance as 
previously applied algorithms without external field based on the paper by Chatelain [11, 
16,21,22] Moreover, it relies partially on the same ideas (see below). It would be 
interesting to compare these algorithms systematically. 

(i) Reweighting spin flip probabilities 

The first idea is to study the response of an infinitely strong external perturbation 
instead of a weak one and to extract the limit h° — > by an appropriate reweighting 
procedure. Reweighting is a standard technique in Monte-Carlo simulations to enhance 
the efficiency (see e.g. [23]). Starting point is the observation that in heat bath dynamics 
as described above a weak perturbation by a positive localised external field hf x <C 1 
changes the probability to get a positively oriented spin -j iew> — _|_i to lowest order as 

i(l + tanh[/3(/ ij + hf xt) )]) = \(l + tanh[/%]) + (^ cosh- 2 [/%]) hf xt) 

+ 0([hf xt) } 2 ). (A.4) 

Therefore, we may equivalently turn up the spin Oj with 100% probability and 
compensate this intervention by weighting the resulting response in the ensemble average 
by the corresponding first-order coefficient. Thus, the reweighting scheme consists of 
the following steps: 

(i) Select random site i at time s. 

(ii) Define a weight Wi = ^/3cosh~ (/3hi), where hi = ^. o~j is the local field caused by 
the nearest neighbours of site % at time s, and store them in the memory. 

(iii) Set o~i := +1 by force. 

(iv) Measure how much this modification increases the magnetisation at site j at some 
later time t > s. 

(v) Multiply the measured response by the weight factor Wi and then perform the 
average over many independent realisations. 
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(ii) Simultaneous tagging of all possible spin flips caused by an external field 

Even with the reweighting procedure described above the measurement of the auto- 
response function is still difficult to perform since the response to a perturbation at a 
single site is quite small and fluctuates strongly. However, in the special case of the Ising 
model with heat bath dynamics it is possible to study the response to a perturbation of 
all spins simultaneously. The reason is that for a given set of random numbers used in 
the simulation a positive external field can only increase the magnetisation in the future, 
turning spins up but never turning spins dowr|§|. Consequently, the cluster of positively 
oriented spins generated in a system without perturbation by an external field is always 
a subset of the actual cluster that would have been generated if the perturbation was 
turned on. Thus, for each site i, where the field may have generated a spin flip, it is 
possible to trace and mark the cluster of all spins in the future that would have been 
turned up by this perturbation. 

The other special property of the Ising model with heat bath dynamics, which is 
exploited by the algorithm, is the circumstance that these clusters do not interact. More 
specifically, if the external field had generated two spin flips, the resulting cluster would 
just be the union of the respective individual clusters. This allows one to tag all spins 
at time s by non-interfering labels and to identify the percolation cluster for each of 
the labels (see Fig. IA1J) . The response of a particular site j is then proportional to the 
number of its labels (re- weighted as described above). This trick accelerates simulations 
by a factor of the order N, where N is the total number of sitesjjjj 

The tagging algorithm can be implemented as follows. In addition to the local 
spin variables o~i we attach to each site % a dynamically generated set Si of labels. For 
example, in C++ such a dynamical set is provided by the class set< . . . > S [N] of the 
standard template library STL. Initially all these sets are empty (Si = 0) and the 
simulation evolves as usual. At time s, however, each site i is tagged by an individual 
label Aj and the weight Wi is recorded as described above. Therefore, at time s each set 
Si has exactly one element, namely, Si = {Aj}. During the subsequent simulation the 
spin variables <7, evolve according to the usual heat bath dynamics as if there was no 
perturbation. However, in each update step the following additional steps are carried 
out 

• When updating spin i, first clear the corresponding set by setting Si — and let 
Ui = (J ■ Sj be the union of all tags carried by the nearest neighbours. 

• For all labels A 6 Ui count the number n\ of nearest neighbours which are oriented 
in negative direction aj = —1 and carry the tag A. Obviously, if an external field 
had flipped the spin corresponding to the label A at time s it would have increased 
the local field hj in the present update by 2n\. 

§ As a consequence, the response function R(t, s,r) is strictly positive, cf. Sect. 2. 
|| The trick would also allow one to study n-point correlation functions efficiently. 
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If hj < h\ and hj + 2n\ > h\ add label A to set Si. 



(c) 



The autoresponse function R(t, s) is then proportional to the weighted sum over all sets 
tagged by their original label, averaged over many independent realisations: 



R(t,s)<x/-J2 W 3 8 ^Si) 



(A.5) 
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Figure Al. Illustration of the algorithm for a small system with 4x4 sites and 
periodic boundary conditions without quenched disorder. The figure shows eight 
snapshots at t = 0,1,..., 7. Initially the system is prepared in a random state. 
Throughout the whole simulation the system evolves as an ordinary unperturbed Ising 
model with heat bath dynamics at T = T c , where black boxes represent positive 
spin variables <7, = 1 while white or grey boxes represent negative spin variables. At 
time t = s = 2 all negatively oriented sites (white boxes) are tagged by individual 
labels a,b,...,i and the corresponding weights are recorded (here w a ~ 0.024482, 
Wb = w c = Wd = w g c± 0.110172, and w e = Wf = Wh — Wi ~ 0.220343). Subsequently 
the simulation tracks how the application of an external local field h > would have 
changed the configuration. For example, if an external field would have turned up the 
spin in the left upper corner at time t = s = 2 (marked by label a) it would have 
changed the configuration at t = 8 in such a way that all sites marked by label a (in 
addition to those marked by black boxes) are positive as well. Note that each site may 
carry several labels, requiring a dynamically generated set structure in the code. At any 
time t > s the autoresponse function is the sum over all sites which carry their original 
label assigned at t = s (indicated by grey boxes), multiplied by the corresponding 
weight factor. In the present example the autoresponse at t = 7 is proportional to w e . 
The measured autoresponse has to be averaged over many independent runs. 
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The following source-code fragment illustrates how the algorithm outlined above can be 
implemented on a computer. The 2d lattice is stored linearly (row by row) with spins 
s [i] and indices % running from to N — 1, where N = L 2 is the number of sites. The 
code requires the existence of a few helper functions which are not listed here, namely: 

• InitLatticeO initialises the lattice in a disordered state. 

• NN(i,d) with d — 0, 1,2,3 computes the index of the four nearest neighbours of 
site i. 

• TruncatedDistance(i, j) computes the geometrical distance between sites % and 
j rounded to an integer. 

• h(i) computes the local field at site i, where the disordered couplings enter. 

The following code fragment performs a single run and computes the response R [r] 
at distance r measured at time t 2 caused by a perturbation at time t±. 



#include <set> 
#include <map> 

const int L=750; 

const int N=L*L; 

const int tl=512; 

const int t2=S*64; 

const int T=l ; 



// system size 

// number of sites 

// time of perturbation 

// time at which the response is recorded 

// temperature (k_B=l) 



int s [N] ; 
double w[N] ; 
double R[L] ; 

set<int> label [N] ; 
set<int>: : iterator p; 



// the lattice ss s[i]=+l,-l 

// weight for each perturbed site 

// the response profile measured at time T 

// set of labels at each site 
// iterator over all labels 



InitLatticeO; // Initialize the lattice 

for (int t=0; t<=t2; t++) { // time loop from t=0 to t=t2 

if (t==tl) for (int n=0; n<N; ++n) { // at t==tl: 

if (s[n]==-l) label [n] . insert (n) ; // assign labels 
w[n] =0 .5/T*pow(cosh(h(n)/T) ,-2 .0) ; // assign weights 
} 
for (int v=0; v<N; ++v) { // perform N random seq. updates 
int i=rndint(N); // select random site 0...N-1 
label [i] . clear () ; // delete its labels 
double z=rnd() ; // draw random number in [0,1] 
if (2*z<l+tanh(h(i)/T)){// if the field is strong enough 
s[n]=l; // turn up the spin 
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else { 

s[n]=-l; // else turn it down 

// Compute the magn. field generated by each label 
// at the neighboring spin-down sites: 

map<int ,double> labelfield; 
map<int ,double> : : iterator k; 

// Loop over all neighboring down spins: 

for (int dir=0; dir<4; dir++) if (S(i,dir)==-1) { 

int j = NN(n,dir) ; 

for (p=label [j] . begin () ; p ! =label [j] . end() ; p++) 
labelfield[*p]+=2*J(i,dir) ; 

> 

// Add all labels that would have created a spin flip: 
for (k=labelfield.begin() ; k! =labelf ield.endO ; ++k) 
if (2*z < l+tanh((h(i) + Ok) .second)/!")) 
label [i] . insert ( (*k) .first) ; 
> > > 
// Compute response profile produced by a single run 
for (int r=0; r<L; ++r) R[r]=0; // clear profile 
for (int i=0; i<N; ++i) 

for (p=label [i] .beginO ; p! =label [i] . end() ; p++) 
R[TruncatedDistance(*p, i)] += w[*p]/N; 



// Perform the average over several runs 
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